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Abstract 

The 25 amino acids long, transmembrane fragment of the Influenza virus M 2 protein 

forms a homotetrameric channel that transports protons across lipid bilayers. It has been 

postulated that high efficiency and selectivity of this process is due to gating by four histidine 

residues that occlude the channel lumen in the closed state. Two mechanisms of gating 

have been postulated. In one mechanism, the proton is “shuttled” through the gate by 

attaching to the S nitrogen atom on the extracellular side of the imidazole ring, followed 

by the release of the proton attached to the c nitrogen atom on the opposite side. In 

the second mechanism, the four histidines move away from each other due to electrostatic 

repulsion upon protonation, thus opening the gate sufficiently that a wire of water molecules 

can penetrate the gate. Then, protons are transported by “hopping” along the wire. In 

this paper, both mechanisms are evaluated in a series of molecular dynamics simulations by 

investigating stability of different protonation states of the channel that are involved in these 

mechanisms. For the shuttle mechanism, these are states with all e protonated histidines, 

one biprotonated residue or one histidine protonated in the 5 position. For the gate opening 

mechanism, this is the state in which all four histidines are biprotonated. In addition, a 

state with two biprotonated histidines is considered. For each system, composed of the 

protein channel embedded in phospholipid bilayer located between two water lamellae, a 

molecular dynamics trajectory of approximately 1.3 ns (after equilibration) was obtained. It 

\ 

is found that the states involved in the shuttle mechanism are stable during the simulations. 
Furthermore, the orientations and dynamics of water molecules near the gate are conducive to 
proton transfers involved in the shuttle. In contract, the fully biprotonated state, implicated 
in the gate opening mechanism, is not stable and the channel looses its structural integrity. 
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If only two histidines are biprotonated the channel deforms but remains intact with the gate 
mostly closed. In summary, the results of this study lend support to the shuttle mechanism 
but not to the gate opening mechanism of proton gating in M 2 . 

A Introduction 

Proton transport across cell membranes is an essential process in cell physiology. This process 

is facilitated and regulated by large, transmembrane (TM) proteins, such as ATP synthases or 

bacteriorhodopsin. The structural complexity of these proteins, however, makes it extremely 

difficult to dissect the molecular mechanisms of protein-mediated proton transport through 

the nonpolar membrane interior. It is thus desirable to have a protein model which is 

small, has a well known structural motif, yet operates with the efficiency and control of 

more complex proteins. This has led to the study of the Influenza-A A7 2 protein - a small, 

homotetrameric, voltage-gated ion channel which transports protons with high efficiency 

and selectivity.[32, 44, 40] Each monomer is built of 97 amino acids and contains a single 

TM domain 19 residues long. Not all residues, however, are essential for transport. Active 

channels have been reconstituted from a synthetic peptide containing only a subset of 25 

residues, including the TM region, with no loss in specificity or efficiency. [13] The sequence 

of amino acids in the peptide is given in Fig 1. The remarkable combination of simplicity and 

efficiency makes M 2 an attractive subject for studies on the mechanism of charge transfer 

\ 

across membranes, and a potential target for re-engineering to create a simple proton pump. 

Although no high resolution structural data for M 2 are available to date, recent NMR 
and CD studies of this protein in phospholipid bilayers strongly suggest that the TM region 
is a-helical, and has the quaternary structure of a four-helix bundle[14, 22]. The center of 
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the bundle forms the transmembrane pore. This model, shown schematically in Fig 1, is 
supported by a recent simulation of M 2 in a water-octane lamella. [45] If the external envi- 
ronment is brought to a pH of approximately 5.3, and a suitable pH gradient exists across the 
membrane, protons are transferred through the channel from the external environment thus 
acidifying the interior of the cell. The efficiency of proton transfer is more than 1000 times 
higher than the efficiencies of transferring other small cations, such as Na + .[11, 7] This, in 
turn, points to the presence of a gating mechanism within the channel which allows protons 
to move through the pore but blocks all larger ions. The titratable group implicated in 
the gating is the intraluminal His-37 residue. [33, 43] In the most widely accepted structural 
models of M 2 , a ring of four histidine residues, one from each of the helices, occludes the 
pore, [41, 33, 45] thus providing an ideal control point inside the channel. The histidines are 
oriented approximately perpendicular to the water-membrane interface with the unproto- 
nated ^-nitrogen atoms of the imidazole rings located on the extracellular side of the gate 
and the protonated 6-nitrogen atoms located on the intracellular side. 

The main aim in this study is to gain insight into the molecular mechanism of proton 
transport and gating in M 2 . To date, two mechanisms have been proposed. [41, 33] In both, 
it has been assumed that protons are translocated along the network of water molecules that 
fills the pore of the channel. This is analogous to the mechanism of proton transfer through 

another channel, gramicidin A. [3] The two mechanisms, however, differ substantially in 

\ 

the proposed way by which protons navigate the gate formed by the histidine residues. In 
one mechanism (further referred to as mechanism I), protons are transferred via a “proton 
shuttle” involving a histidine residue. This is schematically shown in Fig. 2a. In the initial 
state, an e-protonated histidine can accept a proton from the external environment on the 
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unprotonated 5-nitrogen atom of the imidazole ring. This results in a biprotonated inter- 
mediate which is positively charged. This intermediate can relax by eliminating either the 
newly acquired 5-proton back into the environment or the e-proton into the opposite side of 
the channel, thus shuttling the proton across the membrane. Then, the system returns to 
the initial state through tautomerization of the imidazole ring. 

In the alternative mechanism, all four histidines acquire an additional proton and, due 
to repulsion between their positive charges, move away from one another thus opening the 
channel. In the open state, a chain of water molecules penetrates the gate to span the 
full length of the pore. Protons can be efficiently transferred along this chain as long as 
the gate remains open. This mechanism, shown schematically in Fig. 2b, will be called 
mechanism II. It has been proposed on the basis of the results of constrained molecular 
dynamics simulations of a model of M 2 , solvated at both ends by water molecules, and 
placed in a dielectric medium parametrized to represent a membrane core. [41] 

In principle, the correct mechanism can be identified through a series of accurate sim- 
ulation studies of the energetics and dynamics of proton transfer through the M 2 channel 
embedded in a hydrated phospholipid bilayer. In practice, however, such an approach might 
not be feasible, and certainly would be extremely demanding. Since both proposed mecha- 
nisms involve breaking and forming chemical bonds, a methodology that combines classical 

statistical mechanics and quantum mechanics would be desirable. Even if the system is 

\ 

treated in a classical approximation, [39] many complications arise; for example, the need 
to explicitly include polarization effects. [37] It is therefore appropriate to initially take a 
simpler approach, based the fact that each of the proposed mechanisms makes specific, but 
different, predictions about the intermediate states of the protein during proton transfer. 
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In this approach, we require that all intermediate states must be stable, and that the ori- 
entation of water molecules inside the channel must be conducive to proton transfer. If a 
candidate mechanism does not fulfill these requirements, it is unlikely that it provides a 
correct explanation of channel action. 

Each of the intermediate states for mechanisms I and II has been investigated by molecu- 
lar dynamics computer simulations. For mechanism I, we studied the system corresponding 
to the initial state consisting of the four His-37 residues in an e protonation (system la). 
This state has already been simulated in a membrane-mimetic octane lamella [45] and, re- 
cently, in a 1-palmitoyl 2-oleoyl sn-glycero 3-phosphocholine (POPC) membrane. [17] Next, 
we considered a mixed state, consisting of three e, and one biprotonatecl histidine (system 
lb). This is an intermediate in which the proton being transferred is located at the gate 
(see also Fig 2a). Finally, we considered a three c, one S configuration (system lc). This 
state results from a successful proton transfer, but precedes tautomerization of the histidine 
residue which restores the initial state. To test mechanism II, we simulated an intermediate 
containing four biprotonated His-37 residues (system 2a). The same system was studied 
previously in octane lamella. [45] In addition, we also simulated a state with two e and two 
biprotonated species (system 2b) to investigate a possibility that partial protonation of the 
histidines may be sufficient to open the gate and establish a stable water wire spanning the 

channel. The simulations performed are summarized in Table 1. 

\ 

After describing the method of calculations, we present the results for system 1, which 
test the mechanism I, followed by the results for system 2, which corresponds to mechanism 
II. Then we discuss implications of our results for the mechanisms of proton transport in 
M 2 . The paper closes with concluding remarks. 
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B Methods 


The protein was represented as an all-atom model terminally blocked by Ac- and -NHMe 
groups on the N- and C-terminal, respectively. The initial structure was a left-handed, coiled- 
coil, a-helical bundle composed of the 19 TM residues (residues Leu-26 through Asp-44) of 
M 2 . This structure was obtained by computer model building. [12] The modeled protein frag- 
ment was subsequently extended to 25 residues per monomer (Ser-22 through Leu-46), which 
coincides with the sequence used by Duff and Ashley, and Kovacs and Cross in their recon- 
stitution and structural studies. [14, 22] This modification was required because preliminary 
simulations of the channel in a water-hexane membrane mimetic showed that the original, 
shorter model was unstable (the monomers diffused away from each other). There were also 
functional reasons to extend the sequence. It contains three sets of charged residues — Asp- 
24, Asp-44, and Arg-45 — which form “charged rings” at the ends of the bundle. This is a 
common motif in transmembrane ion channels that may be important for filtering ions based 
on their charge and orienting the bundle in the membrane. [9] The rings are schematically 
shown in Fig. 1. All residues forming these rings were simulated in the ionized form. Ion- 
ization states of these residues, however, have not been yet determined experimentally and 
may be different than those assumed here. Recent theoretical calculations of the ionizable 
side chains in the channel-forming helices of the nicotinic acetylcholine receptor have led to 
the suggestion that thd residues comprising the charged rings on the extracellular side of 
the channel are protonated at neutral pH, while the intracellular and intermediate rings are 
fully ionized. [1] It is possible that this is also the case for Influenza M 2 protein. It should 
be, however, kept in mind that the viral channel lacks a significant extramembrane domain 
which may influence ionization states of the nearby residues in the acetylcholine receptor. 
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The protein was embedded in a bilayer built of 90 dimyristoylphosphatidylcholine 
(DlVIPC) molecules and placed between two lamellae of water, each containing 2510 
molecules. An appropriate number of sodium ions was added to maintain system neutrality. 
DMPC is a particularly suitable membrane model since the same phospholipids were used in 
experimental studies on M 2 . [22] The dimensions of the simulation box were 65.0 x 65.0 A 2 in 
the x,y-directions, parallel to the interface, and 200.0 A in the ^-direction, perpendicular to 
the interface. This yields an “open” geometry of the system in which water vapor exists in 
equilibrium with the water lamellae. To eliminate undesired edge effects, periodic boundary 
conditions were applied in the three spatial directions. 

Water molecules were represented by the TIP3P model. [19] The DMPC membrane model 
was developed and tested by Berkowitz et. al.[15] In this model, head groups were represented 
in full atomic detail while the CH n groups in the hydrocarbon chains were treated as united 
atoms, the total mass of which was located on the carbon atoms. To describe the protein, 
the AMBER force field of Cornell et a/.[8] was used. The Lennard-Jones parameters between 
different components of the system were obtained from the standard Lorentz-Bertholot com- 
bination rules. 

For each system listed in Table 1, molecular dynamics simulations were performed at the 
physiological temperature of 310 K in the NVT ensemble for approximately 1.3 ns after equi- 
libration. This temperature is 13 K above the phase transition temperature of DMPC[42], 

\ 

which insures that the bilayer is in the liquid crystal state. Equations of motion were inte- 
grated with a timestep of 2 fs. The AMBER 4.1 simulation package[30] with long range elec- 
trostatic interactions included via the particle mesh Ewald (PME) algorithm[l6] was used. 
Including the long range correction was necessary because the phospholipid head groups 
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carry large atomic charges which interact significantly over distances extending beyond the 
primary simulation cell. 

Before starting the production simulations the system was equilibrated. This was done 
in several steps. First, close interatomic contacts resulting in large van der VVaals repulsive 
energies, found in certain parts of the initial model of M 2 , were relieved by performing several 
short simulations (approximately 20 ps each) of the protein in vacuo. To insure that the 
overall geometry of the bundle did not significantly change during relaxation, the positions of 
non-overlapping residues were kept fixed and the atomic velocities were frequently rescaled. 

Once close contacts between atoms were removed, the protein was inserted into an equi- 
librated configuration of a water-DMPC membrane system. Water and membrane molecules 
that overlapped with the protein were removed and Na + ions were added to keep the system 
neutral. Starting from this configuration, the system was simulated for 500 ps at constant 
pressure with independent scaling of the x,y and 2 dimensions of the box. This allowed 
the water-membrane system to relax around the bundle, which was held in place to keep 
the monomers from diffusing away as the membrane equilibrated around them. Once the 
initial equilibration was achieved, the constraints on the bundle were released and the system 
was simulated for an additional 800 ps. Simulations of the same length were performed to 
equilibrate other systems (in different protonation states of His-37 residues). 

C Results 

C.l Stability of system 1 

The stability of the whole system and, specifically, the protein can be characterized by 
following time evolution of several aggregate, structural indices. The packing arrangement 
of the bundle can be conveniently described by four types of parameters: the average number 
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of residues per turn n,, the average rise per residue d l (i=l,2,3,4) of the individual helices, 
the interhelical separation / 2 ,j mm , and the crossing angle between those members of the 
bundle that are in contact. [ 6 ] In the initial model, the values of r/;, and n, for each helix 
are 1.51 A and 3.86, respectively. By comparison, the corresponding values for an ideal 
a-helical domain are 1.5 A and 3.6, [9]. In real systems, d l varies from 1.4 to 2.0 A, and 
n l varies from 3.2 to 3.9. [ 6 ] The pairs of helices which are in contact in the model of M 2 
are (0,1), ( 0 , 2 ), (1,3) and (2,3). In this model, all of them are separated by an inter-helical 
distance Rij m,n of 9.5 A. The crossing angle for all four pairs of helices is 17°. In other 
proteins these parameters can vary considerably. The values of parameters for the initial 
model are given here as a reference rather than as targets. Because the model is only an 
approximate structural solution, consistent with certain experimental findings rather than 
the precise description of the M 2 structure, deviation from the initial values are expected. 

The time evolution of the four structural properties are shown in Fig. 3. All three 
systems exhibit very similar structural properties, therefore we present the data for system 
la only. The average values for n t , </,, Rij m,n , and f hj are approximately 3.9, 1.5 A, 9.5 A, 
and 38°, respectively. The only quantity that markedly deviates from that in the model is 
the crossing angle. The flux of water in the channel, the degree of solvation of the protein 
ends and the interaction with the membrane are all likely to influence this parameter. It 

is unclear, however, which is the primary factor affecting these changes. The main point 

\ 

emerging from the data shown in Fig. 3 is that the structural features of the bundle are 
stable over the 1 ns. simulation. No systematic drift that would indicate disruption of the 
channel is observed. Substantial fluctuations occur only in the values for n 1 , and n 3 for 
system la. Visual inspection of the protein in this system reveals that helices 1 and 3 are 
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distorted slightly, which yields an overall bend along their long axes. This causes difficulties 
in determining structural parameters because they depend on these axes which, in bent 
structures, can be defined only approximately. 

The conclusion about stability of the three protein systems is reinforced by the data on 
helicity of the monomers as a function of time (not shown). Calculations of this quantity 
using the DSSP program[20] reveal that each of the monomers remain o-helical throughout 
the simulations. The ends of the monomers fray slightly, due to solvation by water, but the 
overall degree of helicity for each monomers in each of the three simulations is greater than 
90% . 

Another measure of the stability of the protein is the calculated root mean square de- 
viation (RMSD) as a function of time. This is shown in Fig. 4. We chose to use the last 
configuration of each trajectory as the reference structure. Proceeding backwards in time 
from the last timepoint, there is an initial adjustment of the structure for all systems, fol- 
lowed by a stable RMSD. No strong drift in RMSD, which would signal a potential structural 
instability, was observed. 

Next we examine the orientation of the helical bundle in the lipid bilayer. In all of the 
simulations the bundle remains, on average, tilted by approximately 30° with respect to the 
bilayer normal. The tilt of individual helices as a function of time in system la is shown in 

Fig. 5. The results for other systems are quite similar. Each helix in the bundle is tilted to 

\ 

a different degree (helix 1 is tilted approximately 10° while helix 3 is tilted 45°), which is 
expected if the bundle as a whole is tilted. In addition, the orientation of each helix in the 
bilayer is stable, i.e., it is not changing markedly during the simulations. 

Tilting of individual a-helices in lipid bilayers has been previously observed in molecular 
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dynamics simulations of polyalanine embedded in a membrane. [42] It was rationalized by 
pointing out that if the hydrophobic core of the membrane is shorter than the length of 
the TM region of the protein, the helix will tilt to maximize its contact with the core. 
Furthermore, our results are consistent with NMR data indicating that the individual helices 
of M 2 exhibit a 30° average tilt relative to the bilayer normal. [22] It should be pointed out, 
however, that the tilt of individual helices does not necessarily imply that the whole bundle 
must be tilted. Since the M 2 geometry is somewhat funnel-shaped, the helices would be tilted 
even if the pore of the channel was aligned with the bilayer normal. In recent simulations 
of M 2 in a membrane mimetic, the bundle did not appear to exhibit significant tilt. [45] 
However, the time evolution of this quantity was not presented but, instead, only a few 
snapshots of the system were shown, making it difficult to ascertain the arrangements of the 
protein during the simulation. 

The discussion of stability, limited so far to the protein, can be further extended to 

include other components of the systems. Good global indices of system stability as a whole 

are density distributions p(Z) of each component measured along the bilayer normal (the 

z-axis). In all cases, they were found not to change as a function of time. The distributions 

for system la, averaged over the full length of the trajectory, are shown in Fig. 6. The 

distributions for systems lb and lc are quite similar and, therefore, are not shown here. In 

each case, the membrane component exhibits well defined head group and hydrocarbon tail 

\ 

regions. The observed structure of the water-membrane interface is consistent with results of 
other simulations of membranes in the absence of integral membrane proteins. [34, 15, 31, 26] 
The only difference is the apparent larger width of the interface in the presence of the protein. 
This is not a result of the deeper penetration of water into the membrane head group region, 
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but reflects the presence of water inside the protein channel. Similarly, the asymmetry of the 
water density profiles on both sides of the membrane is associated with unequal distributions 
of water molecules along the pore of the protein. 

In each system, the protein is symmetrically distributed in the bilayer with its extracel- 
lular N-terminus on the right (z > 0) and intracellular C-terminus on the left (z < 0). The 
overall length of the protein along z is approximately 45 A, and the ends of the protein are 
well solvated by water. By comparing distributions in the two panels of Fig. 6, it can be 
found that charged residues of the protein (Asp-24, Asp-44, Arg-45) are co-located with the 
zwitterionic membrane head groups, which form a highly polar environment. The resulting 
strong electrostatic interactions between these residues and the membrane may help anchor 
the protein in the bilayer and limit the diffusion of the monomers to the plane of the mem- 
brane. In contrast to the charged residues, the gate-forming histidines are located in the 
hydrocarbon core of the membrane, slightly shifted from the center of the bilayer towards 
the intracellular side. The distributions of both charged residues and the histidines do not 
change with time, once again indicating that the system is stable during the simulations. 

C.2 Evaluation of gating predicted by mechanism I 

As already mentioned, it is not only required that all intermediate states must be stable but 
also that the orientation of water molecules inside the channel must be conducive to proton 
transfer. Recall that in mechanism I, a proton reaching the gate becomes attached to a 
His-37 residue on the extracellular side in the 5 position and, subsequently, the e proton is 
released on the intracellular side of the gate. The cycle is completed through tautomerization 
of the histidine residue. Thus, this mechanism implies that the gate remains closed to water 
at all stages of the cycle. Whether this is the case for states la, lb and lc can be examined 
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by calculating the water density profiles in the channel. Since the channel is tilted in the 
bilayer, the channel axis £ is a better coordinate to represent these profiles than the bilayer 
normal. The positions of water molecules are measured relative to the centroid of the His- 
37 residues. As can be seen from the distributions shown in Fig. 7, the channel is nearly 
completely occluded near the gate in all three states. In general, the average number of 
water molecules is small in most of the transmembrane portion of the channel but increases 
markedly toward the ends, as the channel widens. 

Although the number of water molecules in the gate region is quite small, it is not 
precisely zero. Thus, it may be suggested that occasionally a water wire forms across the 
gate that is sufficiently long-lived to transfer proton across the occlusion zone, providing a 
transport mechanism alternative to the “proton shuttle” mediated by a histidine residue. 
This hypothesis can be tested by calculating the persistence time of water clusters spanning 
the gate. 

The calculations required several steps. For each timestep, the water molecules, whose 
oxygen atoms were within 15 A from the channel axis, were selected. Among them were all 
water molecules inside the channel but also some that were located near its edges. Next, 
for each selected water molecule, its neighbors within an oxygen-oxygen cutoff distance of 
3.5 A were identified. This allowed for the construction of an adjacency table whose i- th, 

j-th entry was 1 if water molecules i and j were neighbors or 0 otherwise. Such table is a 

\ 

representation of a directed graph[ 23], and thus water molecules could be assigned to clusters 
by a process known as graph traversal. We implemented a depth-first algorithm to traverse 
the graph. [23] If all water molecules within a 6.0 A cutoff from the His-37 centroid were found 
to belong to one cluster, then the gate was considered to be open. If the water molecules 
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were split between two or more clusters then the gate was considered to be closed. Once 
a water wire across the gate was identified, its persistence time was calculated from the 
number of consecutive configurations (timesteps) in which it was observed. 

The number of events where a water wire formed across the gate, and persisted for at 
least time t ps are shown in Fig. 8. For all systems, these values decay to almost zero in 
less than 4 ps, indicating that there were no water wires which lasted for longer than this. 
In fact, the longest-lived water wire observed in the simulations lasted for slightly less than 
7 ps. Both the scarcity, and the short lifetime of water wires across the gate suggest that 
such wires do not provide an efficient mechanism for proton transport if only one histidine 
undergoes protonation. 

Water molecules near the gate can be further characterized by calculating the radial 
distribution functions, P(r ZJ ), between the nitrogen atoms in the His-37 imidazole ring, and 
water molecules in the channel within a 10 A cutoff. Four different distributions can be 
defined for each histidine residue — two N-0^ a * and two N-H^ a< correlations (using either 
the e or S nitrogen atom of a His-37 as the origin). Note that system la contains four His-37 
residue of one type only. They are in an c protonation state and will be further abbreviated 
as HIE. In system lb, there are three such residues, and one double protonated histidine 
(HIP) carrying a net positive charge. In system lc, there are three HIE and one & protonated 

histidine (HID) residue. The protonation states for each system are listed in Table 1. 

\ 

All four radial distribution functions for the three types of histidine residues are shown 
in Fig. 9. All P(r;vo) exhibit a peak at a typical N-0 hydrogen bonding distance of 2.86 A, 
which indicates that water molecules in contact with the His-37 residues form a well defined 
solvation shell around the imidazole nitrogen atoms. In system la, the peak in P(r^o) is 
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small, indicating that the average number of water molecules residing near each 5-nitrogen 
atom of HIE is relatively small. In fact, the corresponding hydration number, obtained by 
integrating P(r ^o) from zero to the first minimum located at 3.5 A, is equal to 1.0. The first 
peak in P(r i v 4 o) is even smaller and the corresponding hydration number is equal to 0.55. 
For both types of nitrogen atoms, there is also a small peak at 4.8 A. The 3 A separation 
between the first and the second peak is close to a typical oxygen-oxygen distance between 
two hydrogen-bonded molecules of water. 

The suggestion that 5-nitrogen atom of HIE forms a hydrogen bond with the neighboring 
water molecules is confirmed by examining P(r J v (J //). This distribution exhibits peaks at 1.96 
and 3.26 A from the nitrogen atom, consistent with an arrangement in which one hydrogen 
atom of the water molecule is approximately aligned with the N-0 axis and the other forms 
an angle of approximately 105°. In contrast, P(r yv e //) exhibits only one peak at 3.5 A which 
indicates that both hydrogen atoms of the water molecule point away from the N c -H group 
of HIE. The distribution functions for HIE in systems lb and lc are similar and, therefore, 
are not shown here. The only notable difference is that the first peak in P(ryv^o) in system 
lb is somewhat smaller than in the other two systems. 

The radial distribution functions around HID in system lc are also similar. The main 
difference is that the patterns of hydration of 5- and c-nitrogen atoms are reversed, compared 

to HIE, which is consistent with the fact that, in HID, the 5-nitrogen atom is protonated 

\ 

while the e-nitrogen atom is not. Additionally, the second peaks in the distribution functions 
are somewhat larger, indicating a better long range ordering of water molecules. 

The distribution functions for HIP in system lb clearly differ from those for the other 
two types of histidine residues. The distributions around S- and c-nitrogen atoms are nearly 
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identical and suggest the existence of N-H- • • 0 hydrogen bonds, with both hydrogen atoms 
of the water molecule pointing away from the histidine. This agrees well with the fact that 
both nitrogen atoms in HIP are protonated. Furthermore, the first peaks in the distribution 
functions are much larger than the peaks in distributions around nitrogen atoms in the 
other two types of histidine residue (HIE and HID). The corresponding hydration numbers 
are nearly identical (2.0, and 1.9 respectively) since both 5 and e positions are approximately 
equivalent. In addition, the second peaks are very well defined. This shows that the presence 
of a charge on the histidine exerts an attracting influence on water molecules, and induces 
an order in water molecules in the channel that extends for at least two hydration shells. 

To confirm our deductions regarding water orientation around the His-37 residues we cal- 
culated distributions of the angle N — O wat — H wat for each type of imidazole nitrogen. These 
distributions are shown in Fig. 10. For all unprotonated nitrogen atoms, the orientational 
distributions exhibit two clear peaks located at angles of 0° and 104°. This means that one 
hydrogen atom in the water molecule points directly at the unprotonated site along the N-0 
axis. It further indicates that such orientation is fairly rigid. The orientational distributions 
of water molecules near the protonated nitrogen atoms are much broader, suggesting that 
these water molecules undergo a wide range of motions. The common characteristic of these 
motions is that hydrogen atoms preferentially point away from the histidine. It appears that 
the most common arrangement of water molecules around HIP is such that its C 2 v vector 
lies along the N-0 axis and points away from the nitrogen atom (the N — O wat — H wal 
angle is approximately equal to 130°). For uncharged HIE and HID residues, there is also a 
significant fraction of water molecules oriented such that one hydrogen atom points directly 
away from the nitrogen atom, forming the N — O wat — H wat angle of 180°. All these results 
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are consistent with the previously discussed radial distribution functions. The orientations 
of water molecules near different types of histidine residues are schematically summarized in 
Fig. 11. 

Proton transfer capabilities depend not only on the correct orientation of water molecules 
in the channel but also on the lifetimes of these orientations. One way to gauge orientational 
dynamics of a water molecule is by calculating the time correlation functions for its dipole 
and H-H vector defined as: 

Cl(t) = (Pi [r(f + r) • f(r)]) (1) 

where Pi refers to the / th order Legendre polynomial, and r is a vector in the molecular 
frame (either the water dipole or H-H vector). Here, C\ for the dipole vector, and C 2 for the 
H-H vector were calculated by averaging over all time origins r, and over all molecules in a 
specified region. Three types of water molecules were considered — bulk water, water at the 
water-membrane interface, and solvation shell water (within 3.5 A of any His-37 imidadzole 
nitrogen atom). The definion of the His-37 solvation shell is based on the first minimum in 
the N-0 radial distance distributions shown in Fig. 9. 

The results for system la are shown in Fig. 12. The correlations for other systems are 

very similar and thus are not shown. This similarity indicates that the protonation state of 

the His-37 residue does not markedly affect the water dynamics. The results further show 

\ 

that there is a substantial difference in the water dynamics between the three regions. In 
particular, water molecules in the solvation shell of the His-37 residues reorient much slower 
than those in the bulk region and, to a somewhat lesser extent, at the interface. This slower 
motion of water near the His-37 residues could aid proton transfer by keeping the geometry 
of the donor and acceptor molecules fixed. 
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C.3 Evaluation of gating in mechanism II 


In this set of simulations we test the effects of multiple protonations of the His-37 residues. 
System 2a consists of the fully protonated state, whereas system 2b has only two of the His-37 
residues protonated. It is expected that the confinement of the charged His-37 residues to a 
small region in the lumen of the channel will be electrostatically unfavorable in the absence 
of some neutralizing entities, such as counterions. This may cause the histidine residues 
forming the gate to “swing open”, as shown schematically in Fig. 2b, leaving sufficient space 
for water molecules to penetrate the gate. 

For a fully protonated state, containing four HIP residues (system 2a), two molecular 
dynamics trajectories, started from different initial conditions, were obtained. In one tra- 
jectory, the gate was initially in the closed state. In the other trajectory, the protonated 
His-37 residues were first moved to the “open” position characterized in the simulations of 
Sansom et. a/. [41] In both cases the geometry of the helical bundle was initially constrained 
to allow for adjustments of the gate structure without destroying integrity of the channel. 
Once the constraints were released, the fate of the system was the same in both simulations 
— the channel disintegrated within 280 to 450 ps with one monomer separating from the 
remaining three. The separation was preceded by increasing deformation of the bundle, and 
is reflected in the RMSD, as shown in Fig. 13. A representative structure of the disrupted 
bundle is shown in Fig' 14. 

The disruptive effects of electrostatic repulsion should be markedly reduced if the gate is 
only partially protonated. This led us to examine a system in which only two His-37 residues 
are biprotonated (system 2b). The helix parameters for this system, shown in Fig. 15, are 
fairly stable with the exception of n 2 (shown in the top panel). This helix appears to have a 
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bend near the histidine residue, but the regions on both sides of the bend remain a-helical. 
Since the helix is not straight, the helical parameters are less accurate. During the latter part 
of the simulation, the helix becomes less kinked, but simultaneously the histidine residues 
move slightly apart, creating a small gap in the gate region. Overall, the bundle flattens, as 
though the cylinder has been compressed, and looses its conelike geometry. These changes 
are reflected in the RMSD, shown in Fig. 13. Despite rather significant shape changes, 
the monomers remain associated, and provide an adequate channel through the membrane. 
The positioning of the bundle in the membrane is not strongly affected by the changes in 
geometry. 

It is expected that multiple protonation of His-37 residues will increase water penetration 
through the gate region. As can be seen in Fig. 16 for system 2b, the gate is only slightly more 
open relative to that seen in systems 1. Furthermore, the lifetimes of water wires spanning 
the gate remain quite short (see Fig. 8). Both radial and orientational distributions of water 
molecules around HIP and HIE residues are similar to the distributions around the same 
residues in system lb. In general, system 2b undergoes large changes in bundle shape with 
little accompanying changes in features important to the conduction of protons across the 
channel. Thus, it appears that the channel is a robust structure capable of accommodating 
significant changes in bundle geometry. Simultaneous protonation of two His-37 residues 

does not appear to result in excessive amounts of water entering the gate. However, since 

\ 

the simulations do not extend beyond 1 ns, it is not possible to ascertain whether the system 
remains intact over longer time periods. 
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D Discussion 


One of the main objectives of this paper was to examine structural stability of different 
protonation states of the M 2 channel involved in mechanisms I and II. All three protonation 
states required for mechanism I were found to be stable while the channel containing four 
biprotonated His-37 residues, crucial to mechanism II, disintegrated during the simulations. 
These results clearly lend support to the “proton shuttle” mechanism (mechanism I). This 
conclusion, however, should be taken with some caution. Each molecular dynamics trajec- 
tory extended only for approximately 1 ns (excluding equilibration). This is a short time 
compared to timescales of various molecular motions that may influence structures of mem- 
brane embedded proteins. Structures stable for a nanosecond may substantially reorganize 
or even disintegrate over longer periods of time. The reverse, however, is rarely true. If 
a protein assembly rapidly falls apart, it is highly unlikely that it will eventually adopt a 
structure similar to the initial one. Disintegration, of course, may be caused by choosing 
initial conditions for the simulations that are very far from equilibrium — if there are large 
repulsive interactions in the starting configuration any system may rapidly undergo large 
deformations. Precisely for this reason we made considerable efforts to relieve strains in 
system 2a before starting unconstrained, production trajectories. 

Our conclusions about stability of the M 2 channel are in agreement with the results of 
Zhong et al. [45] who studied the same system in a membrane-mimetic lamella of octane. They 
used a different procedure to construct a structure similar to system 2a and also found that 
the helix bundle was unstable. Furthermore, they also observed that the truncated peptide, 
composed of only 19 residues, did not form a proper channel. Both results are at variance 
with simulations reported by Sansom et a/. [41, 17] They found that a system containing four 
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biprotonated His-37 residues was stable and water molecules were able to penetrate the gate. 
In their simulations, however, the membrane was not represented explicitly but modeled as 
a continuum dielectric medium and several artificial restraints were applied to the protein 
bundle. [41] Recently, they reported a 1 ns simulation of the unprotonated, truncated form of 
the channel in the explicit membrane[17] and found this structure to be stable. Since their 
report is very brief, it is not possible at present to determine why their results may differ 
from the results obtained in this study and in the previous work of Zhong et al.. 

We further note that only a few discrete states of the system were considered in this 
work. The actual proton transport also involves a large ensemble of states in which a charge 
resides on a water molecule in the channel. For mechanism I, it is unlikely that these states 
are less stable than system lb, since in this system the charge is located in the most crowded 
part of the channel — the gate. The situation is different in mechanism II. As the proton 
transported along the water wire reaches the gate it brings an additional, fifth positively 
charged species to this region and, thus, the system should become even less stable. 

Another objective of the simulations was to investigate the arrangement of water 
molecules inside the channel and their capability to form hydrogen bonds that assist proton 
translocation. Transport of protons across membranes along networks of hydrogen bonds 
involving water molecules and/or protein side chains is quite common in biological systems. 

It has been postulated to play a role in the functioning of several proteins essential for cellu- 

\ 

lar bioenergetics, such as F^o ATPases, the photosynthetic reaction center, cytochrome / 
and bacteriorhodopsin. [3, 5, 27, 21, 28, 35] The same mechanism has also been implicated 
in proton transport through a simple, transmembrane protein channel, gramicidin A, [3] and 
involves a single file. of water molecules filling the pore. The concept of hydrogen bonded 
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wires is very attractive because it simply explains how protons can translocate through mem- 
branes with high efficiency. The process only requires small displacements of several protons 
between consecutive hydrogen bonded donors and acceptors along the chain, instead of slow 
diffusion of charged species. This idea is essentially the same as in the Grotthuss mechanism 
of proton conductivity in ice. [2] An important characteristic of proton transport along wires 
is that the barrier to the transfer of proton from the donor to the acceptor of a hydrogen 
bond is quite low. In fact, simulations of a model water wire hydrated at both ends showed 
that the free energy of proton transfer is almost constant everywhere along the wire, which 
makes proton transfer nearly barrierless and rapid (on a picosecond time scale). [37] 

The involvement of histidine residues in proton transport via the shuttle mechanism, 
proposed in mechanism I, is also well known in biology. It has been observed in carbonic 
anhydrase, [38, 29] and postulated for respiratory heme-copper oxydases. [28] The ability of 
histidine to tautomerize or change its protonation state in response to the molecular environ- 
ment is well documented for a variety of enzymes. [38, 29, 4, 18] As in proton wire mechanism, 
the barrier to proton transfer via the shuttle mechanism appears to be low. For human car- 
bonic anhydrase III it was estimated at 1.3 kcal/mol. [38] Quantum and statistical mechanical 
studies on model systems, such as NH+-imidazole-NH 3 [25] and HC0O~-imidazole-H20 [24] 
provide a detailed picture of double proton transfer mediated by an imidazole ring. It has 

been shown that concerted mechanism, in which a proton is accepted at the Nj position 

\ 

and another proton is simultaneously donated by the N e atom, is unlikely. Instead, proton 
shuttle proceeds through a double protonated state of the imidazole ring. If both the proton 
donor and proton acceptor are properly aligned with the ring, transport is very fast; in a 
model system the characteristic time for the proton transfer reaction was found to be of the 


22 



order of 1-2 ps. [25] 

These general considerations of proton transport across membranes form the basis for 
discussing implications of our simulations for the mechanism of M 2 action. In both mech- 
anisms considered here, it is assumed that proton translocation proceeds, at least in part, 
through a water wire. Yet, the simulations demonstrate (data not shown) that there is no 
obvious, rigid, hydrogen-bonded chain of water molecules spanning the channel. This does 
not contradict the water wire hypothesis, especially since no charged water species were 
present in the simulated systems, but it points to the highly dynamic nature of the wire. In 
addition, in contrast to gramicidin, water molecules are not arranged in a single file through- 
out the channel. Since the channel has a funnel-like shape, many water molecules fill the 
ends of the pore. In these regions, the orientations of water molecules are influenced by the 
electrostatic field created by highly polar membrane head groups, such that dipole moments 
of water molecules point, on average, towards the middle of the bilayer. This orientation is 
not ideal for a water wire. In fact, the same effect of the electrostatic field will also influ- 
ence orientations of water molecules in gramicidin, but this was not taken into account in 
previous simulations. [36] Considering these arguments, it is possible that hydronium ions 
diffuse through the mouth of the M 2 channel and protons are transported through a proton 
wire/proton shuttle mechanism only through the central part of the pore. 

We further note that water molecules near unprotonated nitrogen atoms of His-37 are 

\ 

oriented such that they form linear 0-H- • • N hydrogen bonds with the imidazole ring. This 
orientation is fairly rigid and relatively long-lived, which is ideal for proton transfer. Simi- 
larly, orientations of water molecules around the N-H groups of the ring are also conducive to 
proton transfer. These arrangements may markedly increase the efficiency of proton translo- 
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cation, compared to other biological systems in which thermodynamic work needs to be ex- 
pended to align all the groups involved in the proton shuttle. For example, proper alignment 
of donor and acceptor groups in human carbonic anhydrase III requires 11 kcal/mol. [38] 

Once a proton is transferred to a ^-nitrogen atom and an approximately symmetric, 
biprotonated HIP state is formed, a question arises: what prevents the backward reaction 
from occurring? Considering that the characteristic time for proton shuttle in model systems 
appears to be quite short, it is possible that a proton does move several times between water 
molecules on the opposite sides of a histidine residue before it is successfully transferred to 
the intracellular side. Since such movements are rapid, they will not have a substantial effect 
on the total efficiency of proton transfer (this is not the rate determining step). Alterna- 
tively, there might be some changes in the position of the histidine residue that decrease the 
probability of the backward reaction but were not captured or identified in the simulations. 
In fact, rotation of His-64 involved in proton shuttle in carbonic anhydrase II was observed 
with the change of pH, [29] but it is not known whether it plays a role in mediating the 
reaction. 

It is natural to assume that the rate determining step in mechanism I is tautomerization 

of His-37 required to bring the system back to its initial state. Unfortunately, there is no 

estimate of the rate of this reaction in M 2 , and no molecular mechanism of tautomerization 

has been proposed. Unquestionably, these issues are central to evaluating the likelihood of 

\ 

mechanism I. We only note that a new cycle of proton transport may start even before the 
previous cycle has been completed. This possibility is implied by the apparent stability of 
system 2b, containing two HIP residues, and by the fact that if a water molecule is located 
near a HfE residue its orientation is unaffected by the presence of the N-H group in another 
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histidine residue on the same side of the gate. Clearly, the ability to process two cycles 
in parallel or nearly in parallel, rather than sequentially, would increase the rate of proton 
transport. 

The nature of the rate determining step clearly distinguishes mechanism I from mech- 
anism II. In the water wire mechanism, this step is most likely the reorientation of water 
molecules in the wire, once a proton has been transferred. For a chain of water molecules 
spanning a membrane, it has been estimated that such a reorientation requires approximately 
500 ps. [26] Note that this step is not of concern in mechanism I; after tautomerization, water 
molecules will spontaneously reorganize to their initial orientations. 

Considering the structural simplicity of M 2 , compared to most other systems transporting 
protons across membranes, understanding the mechanism of its action may be of considerable 
significance for constructing functionalized vesicles. These structures may be very useful in 
medicine and pharmacology, in particular for controlled drug delivery, and may provide 
important models in studies of the origin of cellular life. In particular, it may be possible 
to incorporate into the channel a chromophore capable of ejecting protons in response to 
illumination [10] and re-engineer the protein such that it will act as a simple proton pump, 
an essential component of a bioenergetics system. 

E Conclusions 

\ 

The simulations presented here are consistent with the “proton shuttle” mechanism of proton 
transfer across the closed gate in M 2 , but not with the mechanism that requires the opening 
of the gate and the formation of a continuous water wire. The initial state and two distinct 
intermediates in the “proton shuttle” mechanism that involve different protonation states 
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of one histidine residue were all found to be stable on the timescales of the simulations. 

Furthermore, the orientation of water molecules near the gate with respect to the nitrogen 

atoms in the neighboring histidine residues is conducive to proton transfer across the gate. 

The possibility that a doubly protonated state of the gate can also participate in proton 

transfer is not excluded. In contrast, protonation of all four His-37 residues, which was 

proposed to lead to gate opening, appears to cause structural disintegration of the channel 

due to large, electrostatic repulsion in the gate region. If this repulsion is reduced by allowing 

only two His-37 to be protonated, the gate does not seem to open sufficiently, at least on the 

timescales of the simulations, to yield enough water penetration for efficient proton transport 

through a continuous water wire mechanism. 

Although these results point to a likely candidate mechanism of proton transfer in M 2 , 

many questions regarding this mechanism remain unanswered. What is the free energy and 

kinetics of proton transfer between water molecules and the histidines in the gate? What 

are pK a values of the 8- and c-nitrogen atoms of the histidines at different stages of the 

transport process? What is the likelihood of forward vs. backward transfer of a proton once 

it reaches the gate? What is the time needed for histidine tautomerization? What is the 

role of other residues in the protein in promoting proton transport? To understand fully the 

mechanism of transport, all these questions need to be answered in subsequent experimental 

and theoretical studies. 

\ 
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Table 1: Summary of system configurations — protonation states of His-3 


System 

e-His 

<$-His 

His + 

la 

4 

0 

0 

ib 

3 

0 

1 

lc 

3 

1 

0 

2a 

0 

0 

4 

2b 

2 

0 

2 
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Figure Captions 

Fig. 1: Cartoon of the M 2 bundle with indicators to locations of key residues. The protein 
sequence is indicated below. 

Fig. 2: Sketch of two proposed mechanisms of channel gating by M 2 . Fig. 2a depicts a 
’’proton shuttle” process (mechanism I in this text), whereby proton transfer is mediated by 
the intra-luminal His-37 residues, and is followed by regeneration of the starting configuration 
through tautomerization of the imidazole ring. Fig. 2b depicts a purely structural model 
(mechanism II) in which the open and closed states of the gate are mediated by steric forces 
resulting from protonation of one or more His-37 residues. 

Fig. 3: Helix parameters as a function of simulation time for system la. For the top two 
plots, the residues/turn and the rise/residue are plotted for each individual helix in the 
bundle. For the lower two plots, the separation between helices, and the crossing angles are 
shown for monomers which are in contact. Individual monomer helices are labelled 0, 1,2, 
and 3. 

Fig. 4: Root mean square deviation of the enitire bundle as a function of simulation time 
for systems la, b, and c (legend shown as inset). The RMSD is calculated with respect to 
the last configuration of a given simulation, thus the curves converge to 0 as the end of the 
simulations are approached. 

Fig. 5: Time evolution of the tilt angle of each helix in system 1. The cosine of the angle 
is plotted, where the angle between each helix axis is with respect to the interface normal 
(z-axis). 

Fig. 6: Z dependent density profiles for components of system 1. Upper panel shows water 
and membrane distributions, and total protein distribution. Lower panel shows distributions 
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of selected protein residues. Legends are shown as insets. 

Fig. 7: Number of water molecules as a function of distance along the bundle axis (£). 
Plotted separately are data for systems la (top), lb (middle) and lc(bottom). 

Fig. 8: Probabilities of finding water wires across the gate which persist for at least t 
ps. The details of how this was calculated is described in the text. For all systems, these 
probabilities decay to almost zero in less than 4 ps. The longest-lived water wire observed 
in the simulations lasted for less than 7 ps. 

Fig. 9: Radial probability distributions between imidazole nitrogen atoms and surrounding 
water molecules for His-37. The origin is placed on the nitrogen atoms. Distributions for 
the 5 N are on the left, whereas those for the t N are on the right. 

Fig 10: Orientational distributions for water in the first solvation shell of the imidazole 
nitrogen atoms (as defined by the first minima in the radial distribution functions shown in 
Fig. 9). The solid line corresponds to the angle between defined by 8N — O wat — H wat . The 
light, dashed line shows the angle defined by using the tN in place of the SN. 

Fig. 11: Sketch showing the deduced orientation of water about the His-37 imidazole ring, 
for the different protonation states. 

Fig. 12: Time correlation functions for reorientation of the water dipole, and H-H bond for 
bulk water, water at the membrane surfaces, and in the solvation shell of the His-37 residues. 
Fig. 13: Root mean square deviation of the enitire bundle as a function of simulation time 
for systems 2a, and b. The rest of the caption is the same as for Fig. 4. 

Fig. 14: Snapshot of the M— 2 bundle for the fully protonated system (system 2a) as it 
disintegrates. The carbon backbone is represented as ribbons, and the His-37 residues are 
shown in CPK. 
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Fig. 15: Helix parameters as a function of simulation time for system 2b. The rest of the 
caption is the same as for Fig. 3. 

Fig. 16: Number of water molecules as a function of distance along the bundle axis (£) for 
system 2b. 
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